﻿#include "TestDoseStatCuda.h"

#include <cuda_runtime.h>

#include <chrono>
#include <string>
#include <vector>

#include "Args.h"
#include "DoseStatCuda.h"
#include "MacrosCuda.h"
#include "MyMacro.h"
#include "cnpy.h"

#if 0
const auto& dsInfo = args()->GetTe_DoseStat();
const auto& fPath = fs::path(dsInfo.filePath).generic_string();
#endif

namespace testdosestatcuda
{

    int tt(const mytesttypes::DoseStatIn& dIn)
    {
        MY_ECHOFUNPATH;

        namespace fs = std::filesystem;

#if 0  // in-place 测试数据
        std::vector<double> voxelVolume = {1.0, 2.0, 3.0, 4.0, 5.0};
        std::vector<double> voxelDose = {10.0, 20.0, 30.0, 40.0, 50.0};
        bool doExactAccumulation = true;
#endif

        // 从npz文件读取数据
        std::vector<double> voxelVolume;
        std::vector<double> voxelDose;
        bool doExactAccumulation = false;

        try
        {
            const auto& fPath = dIn.filePath;

            // 对应的 "prescription_type" : "D_95"
            cnpy::npz_t npz_datas = cnpy::npz_load(fPath);
            for (const auto& [key, value] : npz_datas)
            {
                MYINFO("array key: {}", key);
            }
            cnpy::NpyArray arr_volumes = npz_datas[dIn.volumeName];  // 数据存储在键 "volumes"下
            auto* volumes_data = arr_volumes.data<double>();
            voxelVolume = std::vector<double>(volumes_data, volumes_data + arr_volumes.num_vals);
            // 读取剂量数据
            cnpy::NpyArray arr_doses = npz_datas[dIn.doseName];  // 假设数据存储在键"doses"下
            auto* doses_data = arr_doses.data<double>();
            voxelDose = std::vector<double>(doses_data, doses_data + arr_doses.num_vals);

            MYINFO("Success: read data from .npz: ");
            MYINFO("  voxelVolume count: {}", voxelVolume.size());
            MYINFO("  voxelDose count: {}", voxelDose.size());

            // 确保两个数组长度相同
            if (voxelVolume.size() != voxelDose.size())
            {
                MYERROR("Error: voxelVolume and voxelDose have different lengths!");
                // 取较短的长度
                size_t min_size = std::min(voxelVolume.size(), voxelDose.size());
                voxelVolume.resize(min_size);
                voxelDose.resize(min_size);
            }
        }
        catch (const std::exception& e)
        {
            MYERROR("read npz file failed: {}", e.what());
            // 使用默认测试数据
            voxelVolume = {1.0, 2.0, 3.0, 4.0, 5.0};
            voxelDose = {10.0, 20.0, 30.0, 40.0, 50.0};
        }

        MYINFO("Testing Cuda Dose Statistics Functions");
        MYINFO("==================================");

        // 初始化CUDA
        CUDA_CALLC(cudaSetDevice(0));
        MYINFO("CUDA initialization successful");

        // 测试accumulate函数
        MYINFO("\n1. Testing Accumulate function...");
        dosestatcuda::DoseStatRet result = dosestatcuda::Accumulate(voxelVolume, voxelDose, doExactAccumulation);

#if 0
        MYINFO("   Binned Dose Values: ");
        for (const auto& dose : result.binnedDoseValues)
        {
            MYINFO("{} ", dose);
        }
        MYINFO("");

        MYINFO("   Accumulated Relative Volumes: ");
        for (const auto& volume : result.accumulatedRelativeVolumes)
        {
            MYINFO("{} ", volume);
        }
        MYINFO("");
#endif
        ///================================================== 测试get_metric函数
        MYINFO("\n2. Testing GetMetric function...");

        double d98 = dosestatcuda::GetMetric(0.98, result.binnedDoseValues, result.accumulatedRelativeVolumes);
        MYINFO("   D_98 = {} Gy", d98);

        double d95 = dosestatcuda::GetMetric(0.95, result.binnedDoseValues, result.accumulatedRelativeVolumes);
        MYINFO("   D_95 = {} Gy", d95);

        double d50 = dosestatcuda::GetMetric(0.5, result.binnedDoseValues, result.accumulatedRelativeVolumes);
        MYINFO("   D_50 = {} Gy", d50);

        /// 测试get_dose_statistics_roi_interp函数
        MYINFO("\n3. Testing GetDoseStatRoiInterp function...");

        /// 测试平均剂量
        double averageDose =
            dosestatcuda::GetDoseStatRoiInterp("average", result.binnedDoseValues, result.accumulatedRelativeVolumes, voxelVolume, voxelDose);
        MYINFO("   Average Dose = {} Gy", averageDose);

        /// 测试D_95
        double d95_alt =
            dosestatcuda::GetDoseStatRoiInterp("D_95", result.binnedDoseValues, result.accumulatedRelativeVolumes, voxelVolume, voxelDose);
        MYINFO("   D_95 (via interp) = {} Gy", d95_alt);

        /// 测试自定义DoseAtVolume
        double doseAt70 = dosestatcuda::GetDoseStatRoiInterp("DoseAtVolume", result.binnedDoseValues, result.accumulatedRelativeVolumes, voxelVolume,
                                                             voxelDose, 0.7);
        MYINFO("   Dose at 70% Volume = {} Gy", doseAt70);

        return 0;
    }

    ///xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 统计运行时间
    static constexpr int PROFILE_LOOPN_ACCUM = 1000;   // 多次运行取平均值
    static constexpr int PROFILE_LOOPN_METRIC = 1000;  //
    static constexpr int PROFILE_LOOPN_INTERP = 1000;  //
    int tt_profile(const mytesttypes::DoseStatIn& dIn)
    {
        MY_ECHOFUNPATH;

        using namespace std::chrono;

        // 从npz文件读取数据
        std::vector<double> voxelVolume;
        std::vector<double> voxelDose;
        bool doExactAccumulation = false;

        try
        {
            const auto& fPath = dIn.filePath;

            // 对应的 "prescription_type" : "D_95"
            cnpy::npz_t npz_datas = cnpy::npz_load(fPath);
            for (const auto& [key, value] : npz_datas)
            {
                MYINFO("array key: {}", key);
            }
            cnpy::NpyArray arr_volumes = npz_datas[dIn.volumeName];  // 数据存储在键 "volumes"下
            auto* volumes_data = arr_volumes.data<double>();
            voxelVolume = std::vector<double>(volumes_data, volumes_data + arr_volumes.num_vals);
            // 读取剂量数据
            cnpy::NpyArray arr_doses = npz_datas[dIn.doseName];  // 假设数据存储在键"doses"下
            auto* doses_data = arr_doses.data<double>();
            voxelDose = std::vector<double>(doses_data, doses_data + arr_doses.num_vals);

            MYINFO("Success: read data from .npz: ");
            MYINFO("  voxelVolume count: {}", voxelVolume.size());
            MYINFO("  voxelDose count: {}", voxelDose.size());

            // 确保两个数组长度相同
            if (voxelVolume.size() != voxelDose.size())
            {
                MYERROR("Error: voxelVolume and voxelDose have different lengths!");
                // 取较短的长度
                size_t min_size = std::min(voxelVolume.size(), voxelDose.size());
                voxelVolume.resize(min_size);
                voxelDose.resize(min_size);
            }
        }
        catch (const std::exception& e)
        {
            MYERROR("read npz file failed: {}", e.what());
            std::exit(-1);
        }

        MYINFO("Testing Cuda Dose Statistics Functions");
        MYINFO("==================================");

        // 初始化CUDA
        CUDA_CALLC(cudaSetDevice(0));
        MYINFO("CUDA initialization successful");

        ///xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx accumulate time profile
        MYINFO("\n1. Testing Accumulate function...");
        double totalTimeAccum = 0.0;
        for (int i = 0; i < PROFILE_LOOPN_ACCUM; i++)
        {
            const auto iStt = high_resolution_clock::now();

            dosestatcuda::DoseStatRet result = dosestatcuda::Accumulate(voxelVolume, voxelDose, doExactAccumulation);

            const auto iEnd = high_resolution_clock::now();
            const auto duration = duration_cast<microseconds>(iEnd - iStt).count();
            totalTimeAccum += duration;
        }
        const double avgTimeAccum = totalTimeAccum / double(PROFILE_LOOPN_ACCUM);
        MYINFO("\tAccumulate: average execution time ({} runs): {:.2f} microseconds", PROFILE_LOOPN_ACCUM, avgTimeAccum);

        ///xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 测试get_metric函数
        MYINFO("\n2. Testing GetMetric function...");
        dosestatcuda::DoseStatRet result = dosestatcuda::Accumulate(voxelVolume, voxelDose, doExactAccumulation);
        double totalMetricTime = 0.0;
        for (int i = 0; i < PROFILE_LOOPN_METRIC; i++)
        {
            const auto start = high_resolution_clock::now();

            const double d98 = dosestatcuda::GetMetric(0.98, result.binnedDoseValues, result.accumulatedRelativeVolumes);
            const double d95 = dosestatcuda::GetMetric(0.95, result.binnedDoseValues, result.accumulatedRelativeVolumes);
            const double d50 = dosestatcuda::GetMetric(0.5, result.binnedDoseValues, result.accumulatedRelativeVolumes);

            const auto end = high_resolution_clock::now();
            const auto duration = duration_cast<microseconds>(end - start).count();
            totalMetricTime += duration;
        }

        const double avgMetricTime = totalMetricTime / PROFILE_LOOPN_METRIC;
        MYINFO("\tGetMetric: average execution time ({} runs): {:.2f} microseconds", PROFILE_LOOPN_METRIC, avgMetricTime);

        ///===================================================== 测试get_dose_statistics_roi_interp函数
        MYINFO("\n3. Testing GetDoseStatRoiInterp function...");
        double totalInterpTime = 0.0;
        double averageDose = 0.0, d95_alt = 0.0, doseAt70 = 0.0;
        for (int i = 0; i < PROFILE_LOOPN_INTERP; i++)
        {
            const auto start = high_resolution_clock::now();

            averageDose =
                dosestatcuda::GetDoseStatRoiInterp("average", result.binnedDoseValues, result.accumulatedRelativeVolumes, voxelVolume, voxelDose);
            d95_alt = dosestatcuda::GetDoseStatRoiInterp("D_95", result.binnedDoseValues, result.accumulatedRelativeVolumes, voxelVolume, voxelDose);
            doseAt70 = dosestatcuda::GetDoseStatRoiInterp("DoseAtVolume", result.binnedDoseValues, result.accumulatedRelativeVolumes, voxelVolume,
                                                          voxelDose, 0.7);

            const auto end = high_resolution_clock::now();
            const auto duration = duration_cast<microseconds>(end - start).count();
            totalInterpTime += duration;
        }

        const double avgInterpTime = totalInterpTime / PROFILE_LOOPN_INTERP;
        MYINFO("   Average Dose = {} Gy", averageDose);
        MYINFO("   D_95 (via interp) = {} Gy", d95_alt);
        MYINFO("   Dose at 70% Volume = {} Gy", doseAt70);
        MYINFO("   Average GetDoseStatRoiInterp execution time ({} runs): {:.2f} microseconds", PROFILE_LOOPN_INTERP, avgInterpTime);

        return 0;
    }

}  // namespace testdosestatcuda